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Abstract 

We show that standard algorithms for anisotropic diffusion based on centered differencing (including the recent 
symmetric algorithm) do not preserve monotonicity. In the context of anisotropic thermal conduction, this can lead 
to the violation of the entropy constraints of the second law of thermodynamics, causing heat to flow from regions of 
lower temperature to higher temperature. In regions of large temperature variations, this can cause the temperature 
to become negative. Test cases to illustrate this for centered asymmetric and symmetric differencing are presented. 
Algorithms based on slope limiters, analogous to those used in second order schemes for hyperbolic equations, are 
proposed to fix these problems. While centered algorithms may be good for many cases, the main advantage of 
limited methods is that they are guaranteed to avoid negative temperature (which can cause numerical instabilities) 
in the presence of large temperature gradients. In particular, limited methods will be useful to simulate hot, dilute 
astrophysical plasmas where conduction is anisotropic and the temperature gradients are enormous, e.g., coUisionless 
shocks and disk-corona interface. 
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1. Introduction 

Anisotropic diffusion, in which the rate of diffusion of some quantity is faster in certain directions than 
others, occurs in many different physical systems and apphcations. Examples include diffusion in geological 
formations [13] , thermal properties of structural materials and crystals [5] , image processing |ll|4|9j , biolog- 
ical systems, and plasma physics. Diffusion Tensor Magnetic Resonance Imaging makes use of anisotropic 
diffusion to distinguish different types of tissue as a medical diagnostic [2]. In plasma physics, the collision 
operator gives rise to anisotropic diffusion in velocity space, as does the quasilinear operator describing the 
interaction of particles with waves [16] . In magnetized plasmas, thermal conduction can be much more rapid 
along the magnetic field line than across it; this will be the main application in mind for this paper. 

Centered finite differencing is commonly used to implement anisotropic thermal conduction in fusion and 
astrophysical plasmas [6|10|14j . Methods based on finite differencing P and higher order finite elements [15] 
are able to simulate highly anisotropic thermal conduction (x||/X-L 10^; where X|j and x± are parallel 
and perpendicular conduction coefficients, respectively) in laboratory plasmas. "Symmetric" differencing 
introduced in [6] is particularly simple and has some desirable properties: perpendicular numerical diffusion 
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is independent of parallel conduction coefficient x\\i perpendicular numerical diffusion is small, and the 
numerical heat flux operator is self adjoint. While in the symmetric method the components of the heat flux 
are located at cell corners, they are located at the cell faces in the "asymmetric" method. The asymmetric 
method has been used to study convection in anisotropically conducting plasmas [10] and in simulations of 
collisionless accretion disks il4l . 




An important fact that has been overlooked is that the methods based on centered differencing can give 
heat fluxes inconsistent with the second law of thermodynamics, i.e., heat can flow from lower to higher 
temperatures. This accentuates temperature extrema and may result in negative temperatures at some 
grid points, causing numerical instabilities as the sound speed becomes imaginary. Also, in image processing 
applications it is required that no new spurious extrema are generated with anisotropic diffusion [11] , making 
centered differencing unviable. 

We show that both the symmetric and asymmetric methods can be modified so that temperature ex- 
trema are not accentuated. The components of the anisotropic heat flux consist of two contributions: the 
normal term and the transverse term (see §2). The normal term for the asymmetric method (like isotropic 
conduction) always gives heat flux from higher to lower temperatures, but the transverse term can be of 
any sign. The transverse term can be "limited" to ensure that temperature extrema are not accentuated. 
We use slope limiters, analogous to those used in second order methods for hyperbolic problems |19|8|, to 
limit the transverse heat fluxes. For the symmetric method, where primary heat fluxes are located at cell 
corners, both normal and transverse terms need to be limited. Limiting based on the entropy-like function 
(s* = —q ■ VT > 0) is also discussed. 

Limiting introduces numerical diffusion in the perpendicular direction, and the desirable property of 
the symmetric method that perpendicular pollution is independent of x\\ longer holds. The ratio of 
perpendicular numerical diffusion and the physical parallel conductivity with a Monotonized Central (MC; 
see [8] for a discussion of slope limiters) limiter is XJ-,num/X|| ^ 10^"^ for a modest number of grid points (~ 
100 in each direction). This clearly is not adequate for simulating laboratory plasmas which require X||/Xi ~ 
10^ because perpendicular numerical diffusion will swamp the true perpendicular diffusion. For laboratory 
plasmas the temperature profile is relatively smooth and the negative temperature problem does not arise, 
so symmetric differencing [3] or higher order finite elements |15| may be adequate. However, astrophysical 
plasmas can have sharp temperature gradients, e.g., the transition region of the sun separating the hot corona 
and the much cooler chromosphere, or the disk-corona interface in accretion flows. In these applications 
centered differencing may lead to negative temperatures giving rise to numerical instabilities. Limiting 
introduces somewhat larger perpendicular numerical diffusion but will ensure that heat flows in the correct 
direction at temperature extrema; hence negative temperatures are avoided. Even a modest anisotropy in 
conduction (x||/x± 10^) should be enough to study the qualitatively new effects of anisotropic conduction 
on dilute astrophysical plasmas lOj, but the positivity of temperature is absolutely essential for numerical 
robustness. 

The paper is organized as follows. In §2 we describe the heat equation with anisotropic conduction and its 
numerical implementation using asymmetric and symmetric centered differencing. In §3 we present simple 
test problems for which centered differencing results in negative temperatures. Limiting as a method to 
avoid unphysical behavior at temperature extrema is introduced in §4 & §5. Slope limiters are discussed in 
§4 and limiting based on the entropy-like condition in §5. Some mathematical properties of limited methods 
are discussed in §6. In §7 we compare different methods and their convergence properties with some test 
problems. We conclude in §8. 

2. Anisotropic thermal conduction 

Thermal conduction in plasmas with the mean free path much larger than the gyroradius is anisotropic 
with respect to the magnetic field lines; heat flows primarily along the field lines with little conduction in 
the perpendicular direction [3]. In such cases, a divergence of anisotropic heat flux is added to the energy 
equation. Thermal conduction can modify the characteristic structure of the magnetohydrodynamic (MHD) 
equations making it difficult to incorporate into upwind methods. However, thermal conduction can be 
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Fig. 1. A staggered grid with scalars Sij (e.g., n, e, and T) located at cell centers. The components of vectors, e.g., b and q, are 
located at cell faces. Note, however, that for the symmetric method the primary heat fluxes are located at the cell corners [6], 
and the face centered flux is obtained by interpolation (see §2.2). 

evolved independently of the MHD equations using operator splitting, as done in lOJ. The equation for the 
evolution of internal energy density due to anisotropic thermal conduction is 



q = -bn{x\\ - x±)V||r - nx±VT, 



(1) 
(2) 



where e is the internal energy per unit volume, q is the heat flux, x\\ s-nd x± are the coefficients of parallel 
and perpendicular conduction with respect to the local field direction (with dimensions i^T~^), n is the 
number density, T = (7 — l)e/n is the temperature with 7 = 5/3 as the ratio of specific heats for an ideal 
gas, b is the unit vector along the field line, and V|| — b ■ V represents the derivative along the magnetic 
field direction. Throughout the paper we use 7 = 2 to avoid factors of 2/3 and 5/3; results of the paper are 
not affected by this choice. 

We consider a staggered grid with the scalars like n, e, and T located at the cell centers and the components 
of vectors, e.g., b and q, located at the cell faces [17], as shown in Figure[T] The face centered components of 
vectors naturally represent the flux of scalars out of a cell. All the methods presented here are conservative 
and fully explicit. It should be possible to take longer time steps with an implicit generalization of the 
schemes discussed in the paper, but the construction of fast implicit schemes for anisotropic conduction is 
non-trivial. 

In two dimensions the internal energy density is updated as follows. 



At 



+ 



Ax Ay 
where the time step At, satisfies the stability condition [12] (ignoring density variations) 
min[Aa;^, Ay^] 



At < 



2(X|| 



(3) 



(4) 



Aa; and Ay are grid sizes in the two directions. The generalization to three dimensions is straightforward. 

The methods we discuss differ in the way heat fluxes are calculated at the faces. In rest of the section 
we discuss the methods based on asymmetric and symmetric centered differencing as discussed in 6]. From 
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Fig. 2. This figure provides a motivation for using a harmonic average for nx- Consider a 1-D case with the tem- 
peratures and nx's as shown in the figure. Given T_i and T\, and the nx's at the faces, we want to calcu- 
late an average nx between cells —1 and 1. Assumption of a constant heat flux gives, q_i/2 = Qi/2 = 9i i-s., 
— (nx)_i/2(Tb — T_i)/Aa; = —{nx)\/2(Tl ~ To)/ Ax = — — T„i)/2Ax. This immediately gives a harmonic mean, 
which is weighted towards the smaller of the two arguments, for the interpolation nx- 

here on x represent parallel conduction coefficient in cases where an explicit perpendicular diffusion is 
not considered (i.e., the only perpendicular diffusion is due to numerical effects). 

2.1. Centered asymmetric scheme 



The heat flux in the x- direction (in 2-D), using the asymmetric method is given by 



9x,j+i/2,j = -nxbx 



' dT —dT 



(5) 



where overline represents the variables interpolated to the face at {i + f /2, j). The variables without an over- 
line are naturally located at the face. The interpolated quantities at the face are given by simple arithmetic 
averaging, 



{by,i,j-l/2 + i>yA+l,j-l/2 + ^y,i,j + l/2 + j + 1/2 )/ 4, (6) 

9T/ay=(r,,,+i + T,;+i,,+i - T,.,_i - T,+i,,_i)/4Ay. (7) 

We use a harmonic mean to interpolate the product of number density and conductivity [7], 
2 f 1 



nx inx)i,] (nx) 



(8) 



this is second order accurate for smooth regions, but nx becomes proportional to the minimum of the two 
nx's on either side of the face when the two differ significantly. Figure [2] gives the motivation for using 
a harmonic average. Physically, using a harmonic average preserves the robust result that the heat flux 
into a region should go to zero as the density in that region goes to zero, as in a thermos bottle using 
a vacuum for insulation. Harmonic averaging is also necessary for the method to be stable with the time 
step in Eq. (|4]). Instead, if we use a simple mean, the stable time step condition becomes severe by a 
factor ^ max[ni+i j, j]/2min[n,;_|_i j, n^j], which can result in an unacceptably small time step for initial 
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conditions with a large density contrast. Physically, this is because the heat capacity is very small in low 
density regions, so even a tiny heat flux into that region causes rapid changes in temperature. Analogous 
expressions can be written for heat flux in other directions. 



2.2. Centered symmetric scheme 

The notion of symmetric differencing was introduced in [6] , where primary heat fluxes are located at the 
cell corners, with 



9x,i+i/2j+i/2 — -nxhx 



— dT —dT 

OX oy 



(9) 



where overline represents the interpolation of variables at the corner given by a simple arithmetic average, 

bx = {bx,i+l/2,j + ^2;,i+l/2j + l)/2, (10) 
K = (^a,ij + l/2 + ^K,i+lj + l/2)/2, (11) 



dT/dx = (r,+ij + - T,,, - T,,,+i)/2Ax, (12) 

dTjdii=[T,,j+i + - T,- , - T,+ij)/2A2/. (13) 

As before (and for the same reasons), a harmonic average is used for the interpolation of nx, 
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(14) 



Analogous expressions can be written for qy^i+i/2.j+i/2- The harmonic average here is different from [5], who 
use an arithmetic average. Ref. 6 is primarily interested in magnetic fusion applications, where density 
variations are usually well resolved (shocks arc usually not important in magnetic fusion) so arithmetic 
averaging will work well. But there might be some magnetic fusion cases, such as instabilities in the edge 
region of a fusion device, where there might be large density variations per grid cell and a harmonic average 
could be useful. All of the test cases in used a uniform density and so will not be affected by the choice 
of arithmetic or harmonic average. 

The heat fluxes located at the cell faces, qx,i+i/2,j a-nd qy.ij+1/2, to be used in Eq. ([3]) are given by an 
arithmetic average, 

Qx,i+l/2,j = {<lx,i+l/2,j + l/2 + '?x,i+l/2j-l/2 )/2, (15) 
<lyA,j+l/2 = {qyS+l/2.j + l/2 + '7i/,i-l/2 j-+l/2)/2. (16) 

As demonstrated in the symmetric heat flux satisfles the self adjointness property (equivalent to s* = 
—q ■ VT > 0) at cell corners and has the desirable property that the perpendicular numerical diffusion 
(X-L,num) is independent of X||/xj_ (see Figure 6 in fB]). But, as we show later, both symmetric and asymmetric 
schemes do not satisfy the crucial local property that heat must flow from higher to lower temperatures, the 
violation of which may result in negative temperature with large temperature gradients. 

The heat flux in the x- direction qx consists of two terms: the normal term qxx = —nxb^dT / dx and the 
transverse term qxy — —nxbxbydT / dy . The asymmetric scheme uses a 2 point stencil to calculate the normal 
gradient and a 6 point stencil to calculate the transverse gradient, as compared to the symmetric method 
that uses a 6 point stencil for both (hence the name symmetric). This makes the symmetric method less 
sensitive to the orientation of coordinate system with respect to the fleld lines. 

A problem with the symmetric method which is immediately apparent is its inability to diffuse away a 
chess-board temperature pattern as dT/dx and dT/dy, located at the cell corners, vanish for this initial 
condition (see Figure [3]). 
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Fig. 3. The symmetric method is unable to diffuse a temperature distributed in a chess-board pattern. The plus 
(+) and minus (— ) symbols denote two unequal temperatures. The average of dT j ^x\^^x|2,i = (^+ ~ T—)/llx and 
dT/dx\i_^^/2,j+l = (T- - to calculate 9T/aa;|i+y2,i+l/2 = 9^ /dx\^^^i2^j + 9T/aa;|i+i/2,j+l vanishes, similarly 

3. Negative temperature with centered differencing 

In this section we present two simple test problems that demonstrate that negative temperatures can arise 
with both asymmetric and symmetric centered differencing. 

3.f. Asymmetric method 

Consider a 2 x 2 grid with a hot zone (T = f 0) in the first quadrant and cold temperature (T — 0.1) in the 
rest, as shown in Figure^ Magnetic field is uniform over the box with bj. — —by — 1/V2. Number density is 
a constant equal to unity. Reflecting boundary conditions are used for temperature. Using the asymmetric 
scheme for heat fluxes out of the grid point (the third quadrant) gives, qx,i-i/2,j = 1y,i,j-i/2 = Oi ^iid 
1x.i+i/2.j — Qy,i.j+i/2 ~ (9.9/8)nx/Aa; (where Ax = Ay is assumed). Thus, heat flows out of the grid point 
(i, j), which is already a temperature minimum. This results in the temperature becoming negative. Figure 
|4] shows the temperature in the third quadrant vs. time for different methods. The asymmetric method 
gives negative temperature (T^.j < 0) for flrst few time steps, which eventually becomes positive. All other 
methods (except the one based on entropy limiting) give positive temperatures at all times for this problem. 
Methods based on limited temperature gradients will be discussed later. This test demonstrates that the 
asymmetric method may not be suitable for problems with large temperature gradients because negative 
temperature results in numerical instabilities. 

3.2. Symmetric method 

The symmetric method does not give negative temperature with the test problem of the previous section. 
In fact, the symmetric method gives the correct result for temperature with no numerical diffusion in the 
perpendicular direction (zero heat ffux out of the grid point {i,j), see Figure |4]). Other methods resulted in 
a temperature increase at {i,j) because of perpendicular numerical diffusion. Here we consider a case where 
the symmetric method gives negative temperature. 

As before, consider a 2 x 2 grid with a hot zone (T = 10) in the flrst quadrant and cold temperature 
(T = 0.1) in the rest; the only difference from the previous test problem is that the magnetic held lines are 
along the x- axis, b^ — 1 and by — (see Figure [5]) • Reflective boundary conditions are used for temperature, 
as before. Since there is no temperature gradient along the fleld lines for the grid point (i,j), we do not 
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Fig. 4. Test problem to show that the asymmetric method can result in negative temperature. Magnetic field lines are along 



the diagonal with bx 



l/\/2. With the asymmetric method heat flows out of the third quadrant which is already 



a temperature mininmm, resulting in a negative temperature Ti j. However due to immerical perpendicular diffusion, at late 
times the temperature becomes positive again. The temperature at is shown for different methods: asymmetric (solid 

line), symmetric (dotted line), asymmetric and symmetric with slope limiters (dashed line; both give the same result), and 
symmetric with entropy limiting (dot dashed line). 
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Fig. 5. Test problem for which the symmetric method gives negative temperature at Magnetic field is along the x- 

direction, bx = i and by = 0. With this initial condition, all heat fluxes into should vanish and the temperature Tjj 

should not evolve. All methods except the symmetric method (asymmetric, and slope and entropy limited methods) give a 
constant temperature Tij = 0.1 at all times. But with the symmetric method, the temperature at becomes negative due 

to the heat flux out of the corner (i — 1/2, j + 1/2). The temperature Tij eventually becomes equal to the initial value of 0.1. 
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expect the temperature there to change. While all other methods give a stationary temperature in time, 
the symmetric method results in a heat flux out of the grid through the corner at (i — 1/2, j + 1/2). 
With the initial condition as shown in Figure [5l the only non- vanishing symmetric heat flux out of (i,j) is, 
Qx,i-i/2.j+i/2 = — (9.9/2)nx/Ax. The only non-vanishing face-centered heat flux entering the box through 
a face is qx,i-i/2,j = — (9.9/4)nx/Aa; < 0; i.e., heat flows out of which is already a temperature 

minimum. This results in the temperature becoming negative at although at late times it becomes 

equal to the initial temperature at This simple test shows that the symmetric method can also give 

negative temperatures (and associated numerical problems) in presence of large temperature gradients. 



4. Slope limited fluxes 



As discussed earlier, the heat flux is composed of two terms: the normal qxx = —nxb'^dT / dx term, 
and the transverse q^y = —nxbxbydT / dy term. For the asymmetric method the discrete form of the term 
qxx — —nxb^dT/dx has the same sign as —dT/dx, and hence guarantees that heat flows from higher to 
lower temperatures. However, q^y = —nxbxbydT / dy can have an arbitrary sign, and can give rise to heat 
flowing in the "wrong" direction. We use slope limiters, analogous to those used for linear reconstruction 
of variables in numerical simulation of hyperbolic systems [19|8j . to "limit" the transverse terms. Both 
asymmetric and symmetric methods can be modified with slope limiters. The slope limited heat fluxes ensure 
that temperature extrema are not accentuated. Thus, unlike the symmetric and asymmetric methods, slope 
limited methods can never give negative temperatures. 



4.1. Limiting the asymmetric method 



Since the normal heat flux term q^x is naturally located at the face, no interpolation for dT/dx is required 
for its evaluation. However, an interpolation at the x- face is required to evaluate dT/dy used in q^y (the 
term with overlines in Eq. [5]). The arithmetic average used in Eq. ([7|) for dT/dy to calculate qxy was found 
to result in heat flowing from lower to higher temperatures (see Figure |4l). To remedy this problem we use 
slope limiters to interpolate temperature gradients in the transverse heat flux term. 

Slope limiters are widely used in numerical simulations of hyperbolic equations (e.g., computational gas 
dynamics; see [1918) ). Given the initial values for variables at grid centers, slope limiters (e.g., minmod, 
van Leer, and Monotonized Central (MC)) are used to calculate the slopes of conservative piecewise linear 
reconstructions in each grid cell. Limiters use the variable values in the nearest grid cells to come up with 
slopes that ensure that no new extrema are created for the conserved variables along the characteristics, 
a property of hyperbolic equations. Similarly, we use slope limiters to interpolate temperature gradients in 
the transverse heat flux term so that unphysical oscillations do not arise at temperature extrema. 

The slope limited asymmetric heat flux in the x- direction is still given by Eq. ([5]), with the same dT/dx as 
in the asymmetric method, but a slope limited interpolation for the transverse temperature gradient dT/dy^ 
given by 



dT 
dy 



LlL 



i+l/2,3 



dT 



dy 



dT 
^ dT 



L 



dy 



i,i+l/2 

dT 



(17) 



where L is a slope limiter like minmod, van Leer, or Monotonized Central (MC) limiter [Sj; e.g., the MC 
limiter is given by 



MC(a, b) = minmod 



2 minmod(a, 6), ■ 



(18) 



where 
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minmod(a, b) — min(a, b) if a, 6 > 0, 
= max(a, b) if a, 6 < 0, 
= if a6 < 0. 

A slope limiter weights the interpolation towards the argument smallest in magnitude, if the arguments differ 
by too much, and returns a zero if the two arguments are of opposite signs. An analogous expression for the 
transverse temperature gradient at the y- face, dT/dx, is used to evaluate the heat flux qy. Interpolation 
similar to the asymmetric method is used for all other variables (Eqs.|6]&[8]). 



4.2. Limiting the symmetric method 



In the symmetric method, primary heat fluxes in both directions are located at the cell corners (see Eq. 
inj. Temperature gradients in both directions have to be interpolated at the corners. Thus, to ensure that 
temperature extrema are not amplified with the symmetric method, both dT/dx and dT/dy need to be 
limited. 

The face-centered qxx,i+i/2,j is calculated by averaging q^x from the adjacent corners, which are given by 
the following slope- limited expressions: 



1xxA+l/2.j + l/2 



Jxx,i+l/2,j-l/2 



-nxblL2 



-nxblL2 



' dT 




dT 




dx 




dx 


i+l/2j + l_ 


' dT 




dT 




dx 


i+l/2,i 


dx 


i+l/2j-l_ 



(19) 



(20) 



where N and S superscripts indicate the north-biased and south-biased heat fluxes. The face centered heat 
flux used in Eq. ^ is qxx.i+i/2,j = (9^.^+1/2.^+1/2 + 1xx.i+i/2.j-i/2) 1"^'^ ^^'^ interpolated quantities 

(indicated with an overline) are the same as in Eq. ([9|). The limiter L2 which is different from standard slope 
limiters is defined as 

L2(a, 6) = (a + 6)/2, if min(aa, a/a) < (a + 6)/2 < max(aa, a/a), 
= min(aa, a/a), if (a + b)/2 < min(aa, a/a), 

=^ niax(aa, a/ a), if (a + 6)/2 > max(aa, a/ a), (21) 

where < a < 1 is a parameter; this reduces to a simple averaging if the temperature is smooth, while 
restricting the interpolated temperature (dT/dx) to not differ too much from dT/dx\i+if2,j (and be of the 
same sign). We choose a = 3/4 for all of the results in this paper. Note that the L2 limiter is not symmetric 
with respect to its arguments. It ensures that qxx.i+i/2,j±i/2 is of the same sign as —dT/dx\i+i/2j; i.e., 
the interpolated normal heat flux is from higher to lower temperatures. This interpolation will be able to 
diffuse the chess board pattern in Figure [3l The transverse temperature gradient is limited in a way similar 
to the asymmetric method; the temperature gradient dT / dy\i^i/2.j is still given by Eq. (fT7|) . Thus if a = 1, 
the limited symmetric method becomes somewhat similar to the limited asymmetric method (though with 
differences in the interpolation of the magnetic field direction and of nx). 



5. Limiting Vifith the entropy-like source function 



If the entropy-like source function, which we define as s* = —q ■ VT (see Appendix [A] to see how this is 
different from the entropy function) is positive everywhere, heat is guaranteed to flow from higher to lower 
temperatures. For the symmetric method, s* evaluated at the cell corners is positive definite, but this need 
not be true for interpolations at the cell faces; thus heat may flow from lower to higher temperatures. An 
entropy-like condition can be applied at all face-pairs to limit the transverse heat flux terms {qxy and qyx), 
such that 
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dT 



dT 



hj + 1/2 



> 0. 



(22) 



The limiter L2 is used to calculate the normal gradients Qxx and Qyy at the faces, as in the slope limited 
symmetric method (see §4.2). The use of L2 ensures that ~'qxx.i+i/2,jdT/dx\i+i/2j > 0, and only the 
transverse terms qxy and qyx need to be reduced to satisfy Eq. ([22]) . That is, if on evaluating s* at all four 
face pairs the entropy-like condition (Eg. l^^ is violated, the transverse terms are reduced to make s* vanish. 
The attractive feature of the entropy limited symmetric method is that it reduces to the symmetric method 
(which has the smallest numerical diffusion of all the methods; see Figure |9]) when Eq. (|22|) is satisfied. The 
hope is that limiting of transverse terms may prevent oscillations with large temperature gradients. 

The problem with entropy limiting, unlike the slope limited methods, is that it does not guarantee that 
numerical oscillations at large temperature gradients will be absent (e.g, see Figures |4] and [7]) . For exam- 
ple, when dT/dx\i^i/2.j — dT/dy\ij^i/2 = 0, Eq. ([^ is satisfied for arbitrary heat fluxes qx,i+i/2,j and 
Qy^i j+i/2- In such a case, transverse heat fluxes Qxy and qyx can cause heat to flow in the "wrong" direc- 
tion, causing unphysical oscillations at temperature extrema. However, this unphysical behavior occurs only 
for a few time steps, after which the oscillations are damped. The result is that the overshoots are not as 
pronounced and quickly decay with time, unlike in the asymmetric and symmetric methods (see Figures |6] 
&[7]). Although temperature extrema can be accentuated by the entropy limited method, early on one can 
choose sufficiently small time steps to ensure that temperature does not become negative; this is equivalent 
to saying that the entropy limited method will not give overshoots at late times (see Figure [7] and Tables 
[T]|4]) . This trick will not work for the centered symmetric and asymmetric methods where temperatures can 
be negative even at late times (see Figure [7]). 

To guarantee that temperature extrema are not amplified, in addition to entropy limiting at all points, 
one can also use slope limiting of transverse temperature gradients at extrema. This results in a method 
that does not amplify the extrema, but is more diffusive compared to just entropy limiting (see Figure [9]). 
Because of the simplicity of slope limited methods and their desirable mathematical properties (discussed 
in the next section), they are preferred over the cumbersome entropy limited methods. 

6. Mathematical properties 

In this section we prove that the slope limited fluxes satisfy the physical requirement that temperature 
extrema are not amplified. Also discussed are global and local properties related to the entropy-like condition 
s* = -q-VT>Q. 

6.1. Behavior at temperature extrema 

Slope limiting of both asymmetric and symmetric methods guarantees that temperature extrema are not 
amplified further, i.e., the maximum temperature does not increase and the minimum temperature does 
not decrease, as required physically. This ensures that the temperature is always positive and numerical 
problems because of imaginary sound speed do not arise. The normal heat flux in the asymmetric method 
{—Wxh^dT I dx) and the L2 limited normal heat flux term in the symmetric method (Eqs . fT9l and [20|) allows 
the heat to flow only from higher to lower temperatures. Thus the terms responsible for unphysical behavior 
at temperature extrema are the transverse heat fluxes qxy and qyx- Slope limiters ensure that the transverse 
heat terms vanish at extrema and heat flows down the temperature gradient at those grid points. 

The operator L{L{a,b), L(c, d)), where L is a slope limiter like minmod, van Leer, or MC, is symmetric 
with respect to all its arguments, and hence can be written as L{a, b, c, d). For the slope limiters considered 
here (minmod, van Leer, and MC), L{a,b, c, d) vanishes unless all four arguments a,b,c,d have the same 
sign. At a local temperature extremum (say at (i, j)), the x- (and y-) face-centered slopes dT/dy\i j^i/2 and 
dT/dy\ij-i/2 (and dT/dx\i^i/2j and dT/dx\i-i/2j) are of opposite signs, or at least one of them is zero. 
This ensures that the slope limited transverse temperature gradients {dT/dy and dT/dx) vanish (from Eq. 

2 2 

ini). Thus, the heat fluxes are qx,i±i/2,] = -Wh dT/dx\.,±i/2,j and qy^^j±i/2 = -Wh dT/dy\^j±i/2 at 
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the temperature extrema, which are always down the temperature gradient. This ensures that temperature 
never decreases (increases) at a temperature minimum (maximum), and negative temperatures are avoided. 



6.2. The entropy-like condition, s* — —q- VT > 



If the number density n remains constant in time, then muhiplying Eq. ([T]) with T and integrating over 
all space gives 



1 



d 



(7 - 1) dt 



uT^dV 



TV • qdV 



q- WTdV 

nx\y\\T\^dV < 0, 



(23) 



assuming that the surface contributions vanish. This analytic constraint implies that volume averaged tem- 
perature fluctuations cannot increase in time. Locally it gives the entropy-like condition s* = —q ■ VT > 0, 
implying that heat always flows from higher to lower temperatures. 

Ref. [6] has shown that the symmetric method satisfies s* = —q ■ VT > at cell corners. The entropy-like 
function s* evaluated at (z + 1/2, j + 1/2) with the symmetric method is 



^i+l/2j + l/2 



dT 

-qx,i+l/2,j + l/2 -g^ 



i+l/2j + l/2 



qyA+1/2,3 + 1/2 



dT 

dy 



(24) 



'i+l/2j + l/2 



Using the symmetric heat fluxes (Eq. [9]) the entropy-like function becomes, 

2 



s* = nxbx 



dT 


2 


-r-2 


dT 


dx 




F nxby 


dy 



^_——dT dT 



■ nx 



dT —dT 
bx-^ + by — 
dx dy 



> 0, 



(25) 



and integration over the whole space implies Eq. ((23| . Although the entropy- like condition is satisfied by the 
symmetric method at grid corners (both locally and globally), this condition is not sufficient to guarantee 
positivity of temperature at cell centers, as we demonstrate in ^3.2[ Also notice that the modification of 
the symmetric method to satisfy the entropy-like condition at face pairs (see ^ does not cure the problem 
of negative temperatures (see Figure S]). Thus, a method which satisfies the entropy- like condition (s* = 
—q - VT > 0) interpolated at some point does not necessarily satisfy it everywhere, implying that unphysical 
oscillations in the presence of large temperature gradients may arise even if the interpolated entropy-like 
condition holds. 

With an appropriate interpolation, the asymmetric method and the slope limited asymmetric methods 
can be modified to satisfy the global entropy-like condition S* = — J q ■ VTdV/V > 0. Consider 



S* = 



-1 



E 



dT 

qx,^+l/2,J 



dT 

+ <lya,j + l/2 -JT- 



(26) 



where and Ny are the number of grid points in each direction. Substituting the form of asymmetric heat 
fluxes, 



S* 



1 



N^Ny 



E 



nxb-x 



dT 



dT 



dx 



dT 



dy 




»,i + l/2 



(27) 



where overlines represent appropriate interpolations. We define 
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G 



dT 



Gya,j + l/2 - V('^X)i,j + l/2^y,ij + l/2 



dx 

dT 
dy 



i+l/2j 



dT 
dy 



Gx,i,j+l/2 = y/nx^x 



dT 
dx 



(28) 
(29) 
(30) 
(31) 



In terms of G's, Eq. ((27|) can be written as 

^* " NJ^ ^ [^1^+1/2 J + J + l/2 

+ GxA+l/2,jGy^i^l/2,j + GxAJ + l/2Gy^ij + i/2\ ■ 

A lower bound on S* is obtained by assuming the cross terms to be negative, i.e., 



(32) 



si,i + l/2 



|G'a;_i_|-i/2j"Gj, i_^l/2,j I 



\GxA,j+l/2GySJ+l/2 I 



(33) 



Now define Gj,,i+i/2j and G^; ^ as follows (the following interpolation is necessary for the proof to 

hold): 



Gx,i,j+l/2 — L{Gx^i+l/2,j: Gx^i-l/2,jj Gx^i+l/2,j+li G xA-l/2,j + l) : 

GyA+l/2,j — L{Gy i .j^l/2, Gy i j^l/2, Gy i^l j^l/2, G ^ + 1 j _ 1 / 2 ) , 



(34) 
(35) 



where L is an arithmetic average (as in centered asymmetric method) or a slope limiter (e.g., minmod, van 
Leer, or MC) which satisfies the property that \L{a, b, c,d)\ < {\a\ + \b\ + \c\ + |d|)/4. Thus, 



5'*> 



NxNy 



Gl 



1 



+ 1/2J ^y,i,j + l/2 ^ 



\Gx,i+l/2,jGy^ij + l/2\ 



+ \Gx,i+l/2,jGyij_i/2 \ + \Gx^i+l/2,jGyi^ij^i/2 \ + \GxA+l/2,jGyi^ij_i/2\ 
+ \Gy,i,j + l/2Gx,i+l/2,j \ + \Gyjj-i-i/2Gx^i-l/2,j \ + \Gyjj-i-i/2Gx^i+l/2,j + l\ 
+ \Gy,i,j + l/2Gx,i-l/2,j + l\] ■ 

Shifting the dummy indices and combining various terms give, 

1 



(36) 



^* ~ NxNy ^ '^^:i + l/2 J + '^y^i J + 1/2 r, [ 



Gx,i+l/2.jGy.i^j + l/2\ 



+ \Gx,i+l/2,jGy^i^j-i/2 \ + \Gx,i+l/2,jGy^i+ij + i/2 \ + \Gx,i+l/2jGyij^ij^i/2\ 
Jj^Yl \{\Gx,i+l/2,] \ - \Gy^iJ + l/2\) + {\Gx,i+l/2,j \ - \GyA,j-l/2\) 



+ {\Gx,i+l/2,j \ - |Gj^,j+lj + l/2|) + {\Gx,i+l/2,j \ - |Gy_i+ij_i/2|) 

>0. 



(37) 



Thus, an appropriate interpolation for the asymmetric and the slope limited asymmetric methods results in 
a scheme that satisfies the global entropy-like condition. A variation of this proof can be used to prove the 
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global entropy condition 5 > by multiplying Eq. ^ with 1/T instead of T (see Appendix although the 
form of interpolation would need to be modified slightly. It is comforting that introducing a limiter to the 
asymmetric method does not break the global entropy-like condition. However, it is important to remember 
that the entropy-like (or entropy) condition satisfied at some point does not guarantee a local heat flow in 
the correct direction; thus it is necessary to use slope limiters at temperature extrema to avoid temperature 
oscillations. 

7. Further tests 

We use test problems discussed in [10] and [15] to compare different methods. The first test problem 
(taken from [10]) initializes a hot patch in circular field lines; ideally the hot patch should diffuse only 
along the field lines, but perpendicular numerical diffusion causes some cross- field thermal conduction. 
Unlike the limited methods, both asymmetric and symmetric methods show temperature oscillations at 
the temperature discontinuity. The second test problem (from |15| ) includes a source term and an explicit 
perpendicular diffusion coefficient (xi- ) • The steady state temperature gives a measure of the perpendicular 
numerical diffusion x±.num- 

7.1. Diffusion of a hot patch in circular magnetic field 

The circular diffusion test problem was proposed in [TO' . A hot patch surrounded by a cooler background 
is initialized in circular field lines; the temperature drops discontinuously across the patch boundary. At late 
times, we expect the temperature to become uniform (and higher) in a ring along the magnetic field lines. 
The computational domain is a [—1,1] x [—1, 1] Cartesian box. The initial temperature distribution is given 

by 

11 13 
T==12 if 0.5<r<0.7 and ^2^^ ^ 12^^' 

= 10 otherwise, (38) 

where r = ^Jjp^+ip and tan — y/x. Fixed circular magnetic field lines centered at the origin are initialized 
and number density (n) is set to unity. Reflective boundary conditions are used for temperature; magnetic 
field and conduction vanishes outside r — 1. The parallel conduction coefficient x = 0.01; there is no explicit 
perpendicular diffusion {xi_ — 0). We evolve the anisotropic conduction equation ^ till time = 200, by 
when we expect the temperature to be almost uniform along the circular ring 0.5 < r < 0.7. In steady 
state (at late times), energy conservation implies that the ring temperature should be 10.1667, while the 
temperature outside the ring should be maintained at 10. 

Figure [H] shows the temperature distribution for different methods at time=200. All methods result in a 
higher temperature in the annulus r G [0.5, 0.7]. The limited schemes show larger perpendicular diffusion (see 
Tables [l]|4] which give errors, minimum and maximum temperatures, and numerical perpendicular diffusion 
at time=200; also see Figure [5]) compared to the symmetric and asymmetric schemes. The perpendicular 
numerical diffusion (xj^.num) scales with the parallel diffusion coefficient x for methods. Notice that for 
Sovinec's test problem (discussed in the next section) where temperature is smooth and an explicit x_l is 
present, perpendicular numerical diffusion for the symmetric method does not increase with increasing x\\- 

The minmod limiter is much more diffusive than van Leer and MC limiters. Both symmetric and asymmet- 
ric methods give a minimum temperature below the initial minimum of 10, even at late times. At late times 
the symmetric method gives a temperature profile full of non-monotonic oscillations (Figure [S]) . Although 
the slope limited fluxes are more diffusive than the symmetric and asymmetric methods, they never show 
undershoots below the minimum temperature. The entropy limited symmetric method gives temperature 
undershoots at early times which are damped quickly, and the minimum temperature is still 10 at late times 
(see Tables [T]|4] & Figure [T]). Entropy limiting combined with a slope limiter at temperature extrema behaves 
similar to the slope limiter based schemes. 
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Method LI error L2 error Loo error Tmax Tmm X-L,num/X|| 



asymmetric 


0. 


.0324 


0. 


.0459 


0. 


.0995 


10, 


.0926 9 


.9744 


0.0077 


asymmetric minmod 


0. 


.0471 


0. 


.0627 


0. 


.1195 


10, 


.0410 


10 


0.0486 


asymmetric MC 


0. 


.0358 


0. 


.0509 


0. 


.1051 


10, 


.0708 


10 


0.0127 


asymmetric van Leer 


0. 


.0426 


0. 


.0574 


0. 


.1194 


10, 


.0519 


10 


0.0238 


symmetric 


0. 


.0114 


0. 


.0252 


0. 


.1425 


10, 


.2190 9 


.9544 


0.00028 


symmetric entropy 


0, 


.0333 


0, 


.0477 


0, 


.0997 


10, 


.0754 


10 


0.0088 


symmetric entropy extrema 


0. 


.0341 


0. 


.0487 


0. 


.1010 


10, 


.0751 


10 


0.0101 


symmetric minmod 


0. 


.0475 


0. 


.0629 


0. 


.1322 


10, 


.0406 


10 


0.0490 


symmetric MC 


0. 


.0289 


0. 


.0453 


0. 


.0872 


10, 


.0888 


10 


0.0072 


symmetric van Leer 


0, 


.0438 


0, 


.0585 


0, 


.1228 


10, 


.0519 


10 


0.0238 



Table 1. Diffusion in circular field lines: 50 X 50 grid 

The errors are based on the assumption that the initial hot patch has diffused to a uniform temperature (T = 10.1667) in the ring 

0.5< r <0.7, and T = 10 outside it. 



Method LI error L2 error Loo error Tmax Ti„i„ X-L,num/X|| 



asymmetric 


0. 


.0256 


0. 


.0372 


0. 


.0962 


10, 


.1240 9 


.9859 


0.0030 


asymmetric minmod 


0. 


.0468 


0. 


.0616 


0. 


.1267 


10, 


.0439 


10 


0.0306 


asymmetric MC 


0. 


.0261 


0. 


.0405 


0. 


.0907 


10, 


.1029 


10 


0.0040 


asymmetric van Leer 


0. 


.0358 


0. 


.0502 


0. 


.1002 


10, 


.0741 


10 


0.0971 


symmetric 


0. 


.0079 


0. 


.0173 


0. 


.1206 


10, 


.2276 9 


.9499 4.1 X 10"^ 


symmetric entropy 


0, 


.0285 


0, 


.0420 


0, 


.0881 


10, 


.0961 


10 


0.0042 


symmetric entropy extrema 


0. 


.0291 


0. 


.0425 


0. 


.0933 


10, 


.0941 


10 


0.0041 


symmetric minmod 


0. 


.0471 


0. 


.0618 


0. 


.1275 


10, 


.0433 


10 


0.0305 


symmetric MC 


0. 


.0123 


0. 


.0252 


0. 


.1133 


10, 


.1406 


10 


0.00084 


symmetric van Leer 


0. 


.0374 


0. 


.0514 


0. 


.1038 


10, 


.0697 


10 


0.0104 



Table 2. Diffusion in circular field lines: 100 x 100 grid 



Method LI error L2 error Loo error Tmax Ti„i„ X-L,num/X|| 



asymmetric 


0. 


.0165 


0. 


.0281 


0. 


.0949 


10, 


.1565 9 


.9878 


0.0012 


asymmetric minmod 


0. 


.0441 


0. 


.0585 


0. 


.1214 


10, 


.0511 


10 


0.0191 


asymmetric MC 


0. 


.0161 


0. 


.0289 


0. 


.0930 


10, 


.1397 


10 


0.0015 


asymmetric van Leer 


0. 


.0264 


0. 


.0407 


0. 


.0928 


10, 


.1006 


10 


0.0035 


symmetric 


0. 


.0052 


0. 


.0132 


0. 


.1125 


10, 


.2216 9 


.9509 1 


.9 X 10"^ 


symmetric entropy 


0, 


.0256 


0, 


.0385 


0, 


.0959 


10, 


.1103 


10 


0.0032 


symmetric entropy extrema 


0. 


.0260 


0. 


.0391 


0. 


.0954 


10, 


.1074 


10 


0.0032 


symmetric minmod 


0. 


.0444 


0. 


.0588 


0. 


.1219 


10, 


.0503 


10 


0.0192 


symmetric MC 


0. 


.0053 


0. 


.0160 


0. 


.0895 


10, 


.1676 


10 


0.0002 


symmetric van Leer 


0. 


.0281 


0. 


.0426 


0. 


.0901 


10, 


.0952 


10 


0.0038 



Table 3. Diffusion in circular field lines: 200 x 200 grid 



Method 


LI error L2 error Loo error 


T 


max 


T 

mm 


X-L,num/ X|| 


asymmetric 


0.0118 


0.0234 


0.0866 


10. 


,1810 9.9898 5.9 x 10"^ 


asymmetric minmod 


0.0399 


0.0539 


0.1120 


10. 


,0629 


10 


0.0115 


asymmetric MC 


0.0102 


0.0230 


0.0894 


10. 


,1708 


10 


6.8 X 10"* 


asymmetric van Leer 


0.0167 


0.0290 


0.1000 


10. 


,1321 


10 


0.0013 


symmetric 


0.0033 


0.0104 


0.1112 


10. 


,2196 9.9504 


8.4 X 10"^ 


symmetric entropy 


0.0252 


0.0384 


0.0969 


10. 


,1144 


10 


0.0027 


symmetric entropy extrema 


0.0253 


0.0383 


0.0958 


10. 


,1135 


10 


0.0026 


symmetric minmod 


0.0401 


0.0541 


0.1124 


10. 


,0622 


10 


0.0116 


symmetric MC 


0.0032 


0.0122 


0.0896 


10. 


,1698 


10 


6.5 X 10"^ 


symmetric van Leer 


0.0182 


0.0307 


0.1026 


10. 


,1260 


10 


0.0013 



Table 4. Diffusion in circular field lines: 400 x 400 grid 




Fig. 6. The temperature at t = 200 for different methods initiaUzed with the ring diffusion problem on a 400 X 400 grid. 
Shown from left to right and top to bottom are the temperatures for: asymmetric, symmetric, asymmetric-MC, symmetric-MC, 
entropy limited symmetric, and minmod methods. Both asymmetric and symmetric methods give temperatures below 10 (the 
initial minimum temperature) at late times. The result with a minmod limiter is very diffusive. The slope limited symmetric 
method is less diffusive than the slope limited asymmetric method. Entropy limited method docs not show non-monotonic 
behavior at late times, but is diffusive compared to the better slope limited methods. 

Strictly speaking, a hot ring surrounded by a cold background is not a steady solution for the ring diffusion 
problem. Temperature in the ring will diffuse in the perpendicular direction (because of perpendicular 
numerical diffusion, although very slowly) until the whole box is at a constant temperature. A rough estimate 
for time averaged perpendicular numerical diffusion (x_L.iium) follows from Eq. p^, 



J{Tf - T.)dV 
Jdt{JV^TdVy 



(39) 



where the space integral is taken over the hot ring 0.5 < r < 0.7, and Ti and Tf are the initial and final 
temperature distributions in the ring. Figure [8] plots the numerical perpendicular diffusion (using Eq.l39p for 
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0.2 




Fig. 7. Minimum temperature over the whole box as a function of time for the ring diffusion test problem: symmetric (dashed 
line), asymmetric (solid line), and entropy limited symmetric (dot dashed line) methods are shown. Initially the temperature 
of the hot patch is 10 and the background is at 0.1. Both asymmetric and symmetric methods result in negative temperature, 
even at late times. The non-monotonic behavior with the entropy limited method is considerably less pronounced; the minimum 
temperature quickly becomes equal to the initial minimum 0.1. The slope limited heat fluxes maintain the minimum temperature 
at 0.1 at all times, as required physically. 

the ring diffusion problem at different resolutions (see Tables [l]|3]) . The estimates for perpendicular diffusion 
agree roughly with the more accurate calculations using Sovinec's test problem described in the next section 
(compare Figures [5] & O ; as with Sovinec's test, the symmetric method is the least diffusive. Table [S] lists 
the convergence rate of X-L,num for the ring diffusion problem evolved with different methods. 

Table 5 

Asymptotic slopes for convergence of x±,num in the ring diffusion test problem 



Method slope 



asymmetric 


1 


066 


asymmetric minmod 





741 


asymmetric MC 


1 


142 


asymmetric van Leer 


1 


479 


symmetric 


1 


181 


symmetric entropy 





220 


symmetric entropy extrema 


282 


symmetric minmod 





735 


symmetric MC 


1 


636 


symmetric van Leer 


1 


587 



To study the very long time behavior of different methods (in particular to check whether the symmetric 
and asymmetric methods give negative temperatures even at very late times) we initialize the same problem 
with the hot patch at 10 and the cooler background at 0.1. Figure [7] shows the minimum temperature 
with time for the symmetric, asymmetric, and entropy limited symmetric methods; slope limited methods 
give the correct result for the minimum temperature (Tmin — 0.1) at all times. With a large temperature 
contrast, both symmetric and asymmetric methods give negative minimum temperature even at late times. 
Such points where temperature becomes negative, when coupled with MHD equations, can give numerical 
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10 ' • • ' 

0.005 0.015 0.025 0.04 

Ax 

Fig. 8. Convergence of X±,num/X|| a-s the number of grid points is increased for the ring diffusion problem. The numerical 
perpendicular diffusion x± num is calculated numerically by measuring the heat diffusing out of the circular ring (Eg. I39I I. The 
different schemes are: asymmetric (A), asymmetric with minmod (v), asymmetric with MC (□), asymmetric with van Leer 
(*), symmetric (+), symmetric with entropy limiting (o), symmetric with entropy and extrema limiting (>), symmetric with 
minmod (-*-), symmetric with MC (x), and symmetric with van Leer limiter (<3). 

instabilities because of an imaginary sound speed. The minimum temperature with the entropy hmited 
symmetric method shows small undershoots at early times which are damped quickly and the minimum 
temperature is equal to the initial minimum (0.1) after time=l. 

7.2. Convergence studies: measuring x_L.num 

We use the steady state test problem described in [15] to measure the perpendicular numerical diffusion 
coefficient, x±,num- The computational domain is a unit square [—0.5,0.5] x [—0.5,0.5], with vanishing 
temperature at the boundaries; number density is set to unity. The source term Q = 27r^ cos(7ra;) cos(7rj/) 
that drives the lowest eigenmode of the temperature distribution is added to Eq. ([T|). The anisotropic 
diffusion equation with a source term possesses a steady state solution. The equation that we evolve is 

|£ = -V • q + Q. (40) 

The magnetic field is derived from the flux function of the form ip oc cos(7ra:) cos(7r?/); this results in 
concentric field lines centered at the origin. The temperature eigenmode driven by the source function Q is 
constant along the field lines. The steady state solution for the temperature is T{x, y) = cos(7rx) cos(7r?/), 
independent of X||. The perpendicular diffusion coefficient x± is chosen to be unity, thus r^^(0,0) gives a 
measure of total perpendicular diffusion: the sum of x± (the explicit perpendicular diffusion) and x±,num 
(the perpendicular numerical diffusion). 

To account for XJ-,num due to the errors in discretization of the parallel diffusion operator, we calculate 
Xiaium = |r^^(0,0) — T;^Q^(0, 0)1, where Tiso(0, 0) is the central temperature calculated by the discretized 
equations at the same resolution in the isotropic limit X|| = X-L- The convention that we use is slightly 
different (and more accurate) than that used ill previous work, x_L,num — |r-i(0,0) - 1|, which effectively 
assumed that isotropic diffusion gives 1130(0,0) = 1 exactly. 

Figure ini shows the perpendicular numerical diffusivity x±,num — |T~^(0, 0) — Tj~q^(0,0)| for X||/x± = 10; 
100 using different methods. The perpendicular diffusion (x±.num) for all methods except the symmetric 
method increases linearly with X|| . This property has been emphasized by [6] to motivate the use of symmetric 
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0.0125 0.0625 0.1 

Ax 



Fig. 9. A measure of perpendicular numerical diffusion Xx.num = \T~^{0,0) — Tr^\ for X\\/x± = 10 (top) and X\\/x± = 100 
(bottom), using different methods. The different schemes are: asymmetric (A), asymmetric with minmod (v), asymmetric with 
MC (□), asymmetric with van Leer (*), symmetric (+), symmetric with entropy limiting (o), symmetric with entropy and 
extrema limiting (>), symmetric with minmod (*), symmetric with MC (x), and symmetric with van Leer limiter (<). The 
numerical diffusion scales with X|| for all methods except the symmetric differencing [6]. 



differencing for fusion applications, which require the perpendicular numerical diffusion to be small for 
X||/x_L ~ 10^. The slope limited methods (with a reasonable resolution) are not suitable for the applications 
which require X||/X-L 2> 10*; this rules out the fusion applications mentioned in [6|15j . However, only the slope 
limited methods give physically appropriate behavior at temperature extrema, thereby avoiding negative 
temperatures in presence of sharp temperature gradients. The error (perpendicular numerical diffusion, 
XJ_,num = |T~^(0, 0) — ^1^0^(0, 0)1) for most methods except the ones based on minmod hmiter, show a 
roughly second order convergence (see Table [5]). 
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Table 6 

Asymptotic slopes for convergence of error Xx.num = \T~^{0,0) — T.~^"'^(0, 0)| 



Method X\\/X± = 10X||/X± = 100 



a,syninietric 


1 809 
l.OUZ 


1 770 


asymmetric minmod 


U.yo 1 4 


u.y4uo 


asymmetric MC 


1.9185 


1.9076 


asymmetric van Leer 


1.706 


1.728 


symmetric 


1.726 


1.762 


symmetric entropy 


2.407 


2.966 


symmetric entropy extrema 


1.949 


1.953 


symmetric minmod 


0.9155 


0.8761 


symmetric MC 


1.896 


1.9049 


symmetric van Leer 


1.6041 


1.6440 



8. Conclusions 

It is shown that simple centered differencing of anisotropic conduction can result in negative temperatures 
in the presence of large temperature gradients. We present simple test problems where asymmetric and 
symmetric methods give heat flowing from lower to higher temperatures, leading to negative temperatures 
at some grid points. Negative temperature results in numerical instabilities, as the sound speed becomes 
imaginary. Numerical schemes based on slope limiters are proposed to solve this problem. 

The methods developed here will be useful in numerical studies of hot, dilute, anisotropic astrophysical 
plasmas |10|14j , where large temperature gradients may be common. Anisotropic conduction can play a cru- 
cial role in determining the global structure of hot, non-radiative accretion flows (e.g., |1I10I14| ). Therefore, 
it will be useful to extend ideal MHD codes used in previous global numerical studies (e.g., [18]) to include 
anisotropic conduction. Slope limiting methods that prevent negative temperature can be particularly help- 
ful in global disk simulations where there are huge temperature gradients that occur between a hot, dilute 
corona and the cold, dense disk. The slope limited method with an MC limiter appears to be the most 
accurate method that does not result in unphysical behavior with large temperature gradients (see Figures 
[n]&[5]). While we have tried a number of possible variations other than the ones described here, there might 
be ways to further improve these algorithms. Future work might explore other combinations of limiters, 
or limiters on combined fluxes instead of limiting the normal and transverse components independently, or 
might explore using higher-order information to reduce the effects of limiters near extrema while preserving 
physical behavior. 

Although the slope and entropy limited methods in the present form are not suitable for fusion applica- 
tions that require accurate resolution of perpendicular diffusion for huge anisotropy {x\\/x± ^ 10^), they are 
appropriate for astrophysical applications with large temperature gradients. A relatively small anisotropy 
of thermal conduction may be sufficient to study the effects of anisotropic thermal conduction [10 . The 
primary advantage of the limited methods is their robustness in presence of large temperature gradients. 
Apart from the simulations of dilute astrophysical plasmas with large temperature gradients (e.g., magne- 
tized collisionless shocks) , monotonicity-preserving methods may find use in diverse fields where anisotropic 
diffusion is important, e.g., image processing, biological transport, and geological systems. 
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Appendix A. Entropy condition for an ideal gas 



The entropy for an ideal gas is given by >S = nVkln{T^^^'^~^^ /n) + const., where n is the number density, 
V the volume, T the temperature, and 7 the ratio of specific heats (= 5/3 for a 3-D mono-atomic gas). The 
change in entropy that results from adding an amount of heat dQ to a uniform gas is 

,^ nVk dT dQ 

do = — = — — . 

7- 1 T T 

We measure temperature in energy units, so fc = 1. The rate of change of entropy of a system where number 
density and temperature can vary in space (density is assumed to be constant in time) is given by 

* . f = -/..?^^ -/.vi^ = /^"-^ ^ 

where wc use an anisotropic heat flux, if = —nxbb ■ VT. and the integral is evaluated over the whole space 
with the boundary contributions assumed to vanish. The local entropy function defined as s = —q- VT/T'^ 
can be integrated to calculate the rate of change of total entropy of the system. 

In the paper we use a related function (the entropy-like function s*) defined as s* = —q ■ VT to limit the 
symmetric methods using face-pairs, and to prove some properties of different anisotropic diffusion schemes. 
The condition —q ■ VT > ensures that heat always flows from higher to lower temperatures. 
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